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ABSTRACT 


The  N,  N2  system  is  studied  to  yield  rate  coefficients  and  cross  sections  for  molecular 
dissociation  under  conditions  important  for  assessing  nonequilibrium  heating  of 
hypersonic  vehicles.  First  principle  calculations  are  used  to  generate  realistic  nuclear 
interaction  potentials,  and  these  are  used  with  accurate  molecular  dynamics  to  yield  the 
fundamental  data  required  for  a  rigorous  treatment  of  nonequilibrium  chemistry. 


1.0  INTRODUCTION 

Space-craft  entering  planetary  atmospheres  must  withstand  significant  heating  as  the  aero  braking 
arising  from  the  interaction  with  the  atmosphere  converts  kinetic  energy  of  the  craft  into  chemical  energy. 
If  the  entry  speed  is  low  enough,  e.g.  less  than  about  11  Km/s,  which  is  the  entry  speed  for  the  space 
shuttle,  the  rates  of  the  chemical  processes  are  sufficiently  fast  that  chemical  equilibrium  is  an  adequate 
model  to  describe  the  role  of  chemistry  in  the  flow  equations  and  to  model  the  heating.  When  chemistry  is 
in  equilibrium,  the  chemical  composition  and  chemical  energy  is  determined  by  a  single  parameter,  the 
local  temperature. 


When  entry  speeds  increase,  the  chemical  time  scales  and  flow  time  scales  become  more  similar,  and 
reliable  predictions  of  heating  are  no  longer  possible  without  explicitly  dealing  with  nonequilibrium 
chemistry.  For  example,  the  Stardust  return  capsule  entered  the  earth's  atmosphere  at  about  12.6  Km/s,  and 
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Figure  1  Stardust  entering  Earth’s  atmosphere 
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parameters.  The  problem  with  this  solution  is  that  the  model  can  contain  many  uncertainties  that  can  be 
extremely  difficult  to  quantify.  One  often  ignored  problem  is  that  the  phenomenological  nature  of  the 
model  causes  the  input  rate  coefficients  to  lose  their  physical  meaning,  thus  it  is  difficult  to  judge  how  one 
relates  data  obtained  from  experiments  designed  to  measure  rate  data  to  the  quantities  required  for  the 
model.  One  is  then  left  with  validation  of  the  models  by  comparing  to  experimental  data  from  laboratory 
experiments  that  only  approximately  mimic  flight  conditions,  or  by  comparing  to  the  very  sparse  flight 
data.  In  either  case,  the  harsh  conditions  involved  make  the  interpretation  and  exact  characterization  of  the 
data  fiendishly  difficult  and  unambiguous  model  improvement  is  impossible. 

Under  the  support  of  the  NASA  Fundamental  Aeronautics  Hypersonic  research  program,  the 
computational  chemists  at  NASA  Ames  have  embarked  on  an  alternate  solution,  namely  the 
characterization  of  nonequilibrium  chemistry  from  first  principles.  In  this  work  we  will  focus  our  attention 
on  the  dissociation  of  N2  by  means  of  collisions  with  N  atoms  and  N2  molecules. 


2.0  CHEMISTRY  BASED  FORMALISM 

The  theoretical  formulation  encompasses  vast  length  scales,  from  macroscopic  vehicle  dimensions  to 
atomistic  dimensions.  The  makes  a  rigorous  complete  formulation  impossible.  Fortunately  there  are 
numerous  physical  approximations  which  are  quite  reliable  that  allow  us  to  decouple  various  length  scales. 
Thus  we  will  take  a  divide  and  conquer  approach. 

We  start  at  the  atomistic  level.  Molecules  and  atoms  move  very  fast  and  the  range  of  inter-molecular 
interactions  is  fairly  short  range,  thus  the  time  between  collisions  is  much  greater  than  the  duration  of  a 
typical  collision.  Thus,  to  a  very  good  approximation,  we  can  start  by  considering  our  universe  to  consist 
of  two  particles  that  undergo  a  collision  in  isolation  from  all  other  particles. 

2.1  N+N 

We  consider  initially  two  nitrogen  atoms.  How  do  we  describe  their  motion?  This  by  itself  is  a  difficult 
task,  principally  because  the  forces  between  particles  is  due  not  only  to  the  nuclei,  but  also  to  the 
electrons.  Thus  we  must  explicitly  treat  the  electrons  and  nuclei  of  the  colliding  particles.  This  leads  to 
high  dimensionality  problems:  for  example,  a  nitrogen  atom  has  seven  electrons,  so  the  treatment  of  two  N 
atoms  is  a  problem  in  3x2x(l+7)=48  dimensions. 

Before  we  go  on  and  describe  how  we  proceed,  it  should  be  noted  that  nuclear  spin  plays  a  non-trivial 
role.  [2]  The  nuclear  spin  leads  to  very  small  corrections  to  the  forces  between  the  particles,  thus  can  be 
ignored  when  determining  the  interactions.  However,  since  a  N  nucleus  has  a  nuclear  spin  1,  it  is  a  boson, 
and  thus  only  spin-electronic-ro-vibrational  wavefunctions  that  are  totally  symmetric  are  present  in  nature. 
Thus  for  the  ground  electronic  state  of  N2,  out  of  the  9  possible  nuclear  spin  couplings,  only  3  appear  with 
odd  rotational  quantum  numbers,  and  only  6  appear  with  even  quantum  numbers.  In  spectra,  this  is 
manifested  by  a  2:1  alternation  in  the  intensities  of  the  different  rotational  levels.  This  complication  on 
the  kinetics  can  be  included  in  a  straightforward,  rigorous  manner,  but  introduces  additional  complications 
into  the  formalism[3]  which  will  be  neglected  in  the  present  work  to  make  the  presentation  more 
transparent. 

We  will  now  speak  of  nuclear  spin  no  more  and  go  on  to  the  treatment  of  the  nuclei  and  electrons.  A 
problem  of  this  sort  is  not  feasible  to  solve  exactly,  thus  we  must  resort  to  further  approximations. 
Fortunately  nature  is  kind,  and  the  ratio  of  N  nuclear  mass  to  an  electron  mass  is  about  2.6xl04,  which 
enables  us  to  obtain  excellent  results  if  we  consider  the  nuclei  to  be  stationary  when  we  treat  the  electrons. 
[2]  For  N2,  this  enables  us  to  reduce  the  problem  to  one  in  3x2=6  dimensions  (the  nuclear  problem),  and 
many  in  3x2x7=42  dimensions  (the  electron  problem). 
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The  electron  problem  remains  a  very  formidable  problem,  but  we  again  have  several  mitigating  factors. 
The  electron  motion  is  governed  by  the  laws  of  quantum  mechanics,  which  simply  put  leads  to  probability 
distribution  functions  for  particles,  and  the  probability  distribution  goes  hand  in  hand  with  an  average 
energy,  with  the  energy  levels  fixed  by  boundary  conditions  and  Planck's  constant.  The  first  mitigating 
factor  is  that  all  electrons  are  indistinguishable,  which  means  that  a  mean  field  approximation  can  provide 
a  reasonable  starting  point  for  obtaining  a  solution.  In  the  mean  field  approximation,  the  motion  of  one 
electron,  a  3  dimensional  problem,  is  found  from  using  the  average  position  of  all  the  remaining  electrons. 
This  leads  to  an  iterative  refinement  of  the  orbitals  describing  the  probability  distribution  of  an  electron's 
position.  This  level  of  calculation  is  quite  feasible  for  very  large  systems  (hundreds  of  atoms).  To  go 
beyond  this  level  of  approximation,  which  is  required  for  the  accuracy  we  need,  we  can  consider  two 
electrons  at  a  time,  with  the  motion  of  the  remaining  electrons  being  averaged  over.  This  is  more 
expensive,  and  is  still  not  accurate  enough.  The  most  reliable  method  used  today  is  to  treat  the  motion  of 
two  electrons  explicitly,  but  allow  the  other  electrons  to  interact  via  functional  restrictions.  Very  efficient 
computer  codes  have  been  written  to  solve  for  this  problem,  but  still  solutions  are  practical  only  up  to  half 
a  dozen  or  so  atoms  in  the  system. 

The  second  mitigating  factor  is  that  the  light  mass  of  the  electrons  means  that  the  spacing  of  energy 
levels  is  very  great.  Thus,  to  a  very  good  degree  of  approximation,  we  need  to  consider  only  the  lowest 
possible  energy  level  of  the  electrons.  Nonetheless,  a  calculation  of  the  electron  energy  level  for  a  single 
nuclear  geometry  requires  the  solution  of  a  non-linear  problem  with  millions  of  unknowns. 

By  splitting  the  nuclear  and  electron  problems,  we  have  enabled  the  solution  of  the  electron  problem 
without  consideration  of  the  nuclear  motion,  but  at  the  expense  of  having  to  solve  repeatedly  the  electron 
problem  for  each  position  the  nuclei  take  during  the  collision.  This  is,  of  course,  another  formidable 
problem,  because  the  nuclei  sample  an  infinite  number  of  positions  during  a  collision.  To  circumvent  this 
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Figure  2:  Calculated  and  fitted  potential  energy  curve  of  N2  and  simple  Morse  model. 

problem,  we  carefully  select  a  set  of  nuclear  geometries  which  we  think  well  characterize  the  set  of 
geometries  encountered,  and  solve  for  the  electronic  energy  at  these  points.  These  points  are  then  fit  to 
some  analytic  function  in  order  to  interpolate/extrapolate  the  energies.  In  Figure  2  we  show  such  a  result 
for  N2.  We  show  the  data,  the  fitting  function,  and  finally  a  Morse  potential.  The  Morse  potential  is  a 
simple  three  parameter  function  which  has  simple  vibrational  energy  levels,  while  the  fit  is  a  14  parameter 
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analytic  function.  The  fit  reproduces  the  data  very  well,  and  the  Morse  function  does  a  good  job  of 
representing  the  data  until  about  half  way  up  to  dissociation  limit,  but  then  starts  to  give  rise  to  large 
deviations.  Since  in  the  present  work,  N2  dissociation  is  of  interest,  it  is  likely  that  the  Morse  potential 
would  not  provide  a  good  basis  for  further  work.  It  is  interesting  to  note  that  while  the  Morse  potential 
appears  to  have  longer  range,  i.e.  it  seems  to  decay  more  slowly  than  the  fit,  the  fit  is  actually  of  longer 
range,  for  it  approaches  the  asymptotic  limit  as  r"6,  where  r  is  the  N2  bond  length,  while  the  Morse 
potential  approaches  the  asymptotic  limit  exponentially.  [2]  This  has  profound  implications  for  the  number 
of  bound  ro-vibrational  levels  of  N2.  Also  note  the  higher  density  of  data  points  in  the  vicinity  of  a  bound 
length  of  4.  This  was  required  to  ensure  the  analytic  representation  behaved  properly  as  the  potential 
evolved  from  its  asymptotic  long  range  behaviour  to  the  short  range  bonding  interaction.  This  effected  the 
number  of  bound  ro-vibational  levels,  and  will  be  discussed  further  below. 

Let  us  now  move  on  to  the  treatment  of  nuclear  motion.  After  the  electrons  have  been  taken  care  of,  we 
have  6  degrees  of  freedom  remaining.  By  transforming  to  a  new  set  of  cartesian  position  vectors,  one  that 
specifies  the  center  of  mass  of  the  N2  molecule,  and  the  other  which  specifies  the  length  and  direction  of 
the  bond,  we  can  simplify  the  problem,  since  the  potential  energy  curve  does  not  depend  on  the  center  of 
mass  cartesians.  This  reduces  the  problem  to  one  in  three  dimensions,  since  the  center  of  mass  motion  is 
easily  solved  for  and  furthermore  plays  no  role  in  this  study.  By  transforming  the  cartesians  specifying  the 
length  and  direction  of  the  bond  into  spherical  polar  coordinates,  we  can  further  simplify  the  problem, 
since  the  potential  energy  curve  only  depends  on  r,  the  length  of  this  vector,  and  partially  separate 
vibration  and  rotational  motion. 

Thus  we  now  have  a  rotational  and  vibrational  problem  to  solve.  How  do  we  do  this?  By  virtue  of  their 
small  mass,  electrons  must  be  treated  by  quantum  mechanics.  As  particles  get  heavier,  gradually  classical 
mechanics  becomes  more  valid.  Now  since  the  mass  of  an  N  atom  is  about  2.6xl04  times  greater  than  an 
electron’s  mass,  it  is  tempting  to  consider  using  classical  mechanics.  This  is  because  classical  mechanics  is 
local  while  quantum  mechanics  is  non-local  and  hence  much  more  expensive.  Thus  let  us  consider  three 
ways  for  treating  nuclear  dynamics:  classical  mechanics,  old  quantum  mechanics,  and  quantum 
mechanics. 


In  a  classical  mechanical  calculation,  one  determines  the  time  evolution  of  a  set  of  coordinates  qu  and 
conjugate  momentum pt.  The  time  evolution  is  governed  by  Hamilton’s  equations[4] 


4>i 


m 


where  the  Hamiltonian  for  a  1-d  vibrator,  for  example,  is 

H  =  +  V  (q) 

2m 

where  m  is  the  reduced  mass  for  the  vibrator,  mN/ 2  for  N2.  The  solution  of  these  equations,  while 
interesting  by  themselves,  do  not  appear  helpful  in  the  situation  of  approximating  quantum  mechanics,  for 
nothing  about  quantization  of  energy  appears.  To  proceed,  we  must  invoke  the  premises  of  so-called  “old” 
quantum  mechanics.  [5]  In  this  procedure,  we  take  the  above  equations  and  transform  to  a  new  set  of 
variables,  called  action-angle  variables.  These  variables  are  special  in  that  the  new  momenta,  the  actions, 
are  constant  in  time,  i.e.  the  new  Hamiltonian  has  no  dependence  on  the  variables  conjugate  to  the  actions, 
the  angle  variables.  Then  the  quantization  rule  is  that  the  only  allowed  values  of  the  actions  are  half 
integral  multiples  of  h,  but  the  particles  still  follow  regular  trajectories:  the  angle  variables  progress  in 
straight  lines.  The  action-angle  variables  can  be  determined  analytically  for  a  few  1-d  problems,  but  must 
be  determined  numerically  for  most  systems  of  interest. 


We  now  compare  this  to  quantum  mechanics.  Now  we  must  solve  a  eigenvalue  eigenvector 
problem: 
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where  the  quantum  mechanical  Hamiltonian  operator  is  built  from  simple  rules  from  the  classical 
Hamiltonian,  IT'J2  gives  the  probability  density,  and  En  is  the  energy.  For  the  1-d  vibrator,  the  Hamiltonian 
operator  is 

h2  #  r//  , 

- T  F  V (<7)- 

2m  dq 

How  do  the  predicted  energies  between  the  old  and  new  quantum  mechanics  compare?  For  a 
diatomic  molecule,  like  N2,  the  agreement  is  quite  good.  For  example,  the  rotational  energy  levels  from 

quantum  mechanics  are  +  ,  where  j  is  the  integral  rotational  quantum  number  and  r  is  the  bond 


2  mr 


(j+\)n2 

length,  while  the  rotational  energy  levels  from  old  quantum  mechanics  are  A - — — .  The  difference 

2  mr 


between  these  two  results  is  only 


r 


which  is  pretty  small  considering  how  large  m  is.  For  vibration, 


8  mr1 

the  situation  is  similar:  for  a  Morse  oscillator  and  Harmonic  oscillator,  the  old  quantum  mechanics  gives 
the  same  answer  as  does  quantum  mechanics,  provided  the  domain  of  the  bond  length  is  -qo  to  +oo,  a  not  to 
restrictive  assumption  for  the  current  value  of  m.  Thus  we  see  that  the  old  quantum  mechanics  works  quite 
well  for  N2. 


Now  consider  what  happens  when  one  couples  vibration  and  rotation.  One  finds  that  the  net  effect 


is  to  add  the  term 


2  mr 


2  to  the  vibrational  potential.  It  should  be  noted  that  the  attraction  of  a  Morse 


potential,  namely  that  action-angle  variables  can  be  determined  analytically,  does  not  hold  when  the 

(j  + 1)2  h2 

— - —2 —  that  is  added  to  the  potential  precludes  analytic  action-angle 


molecule  rotates.  The  term 
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variables. 


Figure  3:  Potential  curve  of  N2  for 7=0,20,40,.. 
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The  effect  of  rotation  on  the  potential  of  N2  is  shown  in  Figure  3.  In  general,  as  j  increases,  the 
potential  minimum  decreases  in  depth  and  gets  pushed  to  larger  r.  An  interesting  situation  occurs  when  the 
potential  minimum  rises  to  the  point  it  is  no  longer  the  global  minimum,  but  simply  a  local  minimum.  In 
this  case,  quantum  mechanically,  there  is  probability  of  finding  the  particle  inside  of  the  centrifugal  barrier 
as  well  as  outside  the  barrier.  Physically,  if  the  particle  started  inside  the  centrifugal  barrier,  the  probability 
of  finding  it  there  will  decay  exponentially  with  time.  These  levels  are  called  quasi-bound  levels, 
predissociated  levels,  or  resonant  levels,  depending  on  the  scientific  field  in  which  they  are  discussed. 
Finally,  eventually  as  j  becomes  larger,  the  potential  will  become  purely  repulsive,  and  no  more  bound  or 
quasi-bound  levels  will  exist. 

In  Figure  4,  we  show  the  range  of  v  and  j  for  which  there  are  bound  or  quasibound  levels.  For  j= 0, 
there  are  bound  vibrational  levels  up  to  60,  whereas  for  v=0,  there  are  bound  rotational  levels  up  to  j  about 
210,  and  quasibound  levels  continuing  up  to  j- 279.  In  our  initial  calculations  on  N2,  we  noticed  a  rather 
odd  bulge  on  a  plot  of  this  type  in  the  vicinity  of  v=50.  This  turned  out  to  be  due  to  additional  bound  states 
coming  from  incorrect  behaviour  of  the  potential  of  Figure  2  in  the  region  around  4  bohr.  When  we  added 
the  additional  points  around  4  bohr  and  improved  the  analytic  representation  in  this  region,  the  bulge  went 
away,  and  we  were  left  with  a  plot  very  like  Figure  4.  The  potential  used  to  generate  the  data  shown  in  this 
figure  is  obtained  not  from  the  potential  of  Figure  2,  but  rather  an  accurate  potential  determined  from 
fitting  spectroscopic  data  for  N2.[6]  The  differences  between  these  potentials  is  not  that  great,  but  certainly 
the  one  based  on  experiment  is  more  accurate,  so  we  will  use  it.  The  number  of  bound  and  quasibound 
levels  for  N2  is  9390  from  old  quantum  mechanics. 


Figure  4:  Rotational-vibrational  levels  of  N2. 

Let  us  consider  the  approximation  of  decoupling  vibration  and  rotation  for  a  moment.  At  the  very 
least,  for  this  approximation  to  be  reasonable,  it  is  necessary  for  the  bound  v,j  states  of  interest  to  fall 
within  a  rectangle  in  a  plot  such  as  Figure  4.  This  occurs  for  v  up  to  about  30  and  j  up  to  about  125.  This 
does  not  include  any  levels  near  dissociation.  When  dissociation  is  important,  the  vj  states  of  interest  fall 
within  a  shape  more  like  a  triangle.  Furthermore,  decoupling  vibration  and  rotation  precludes  the  existence 
of  predissociated  levels.  However,  these  levels  are  expected  to  play  a  very  important  role  in  dissociation. 
Thus  for  dissociation,  the  decoupling  of  vibration  and  rotation  can  not  be  a  reasonable  assumption. 
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2.2  N+N2 

We  now  add  another  particle  to  our  mix,  e.g.  we  consider  our  universe  to  consist  of  N  and  N2.  The 
solution  of  the  electronic  structure  problem  proceeds  as  for  N2,  except  now  things  are  more  complex 
because  of  the  larger  number  of  electrons  and  because  the  potential  energy  hypersurface  now  depends  on 
three  coordinates,  rather  than  one.  One  interesting  complication  is  that  since  the  ground  electronic  state  of 
the  N  atom  is  4S,  the  N3  potential  energy  hypersurface  has  between  3  and  9  open  shell  orbitals,  depending 
on  whether  the  system  looks  more  like  N+N2,  or  3N.  The  large  number  of  open  shell  orbitals  puts 
considerable  stresses  on  the  software  we  use  to  compute  the  potential  energy  hypersurface.  [7]  In  contrast 
to  the  situation  with  N2,  where  we  carried  out  electronic  structure  calculates  and  tens  of  geometries,  for  N3 
we  carried  out  calculations  for  hundreds  of  geometries.  An  analytic  representation  of  these  calculations 
was  then  carefully  formed  and  used  for  all  subsequent  calculations. 

It  is  informative  to  consider  the  general  topology  of  the  potential  energy  hypersurface.  The  N2 
potential  energy  curve  in  Figure  2  has  a  single  minimum  at  r  about  2,  and  then  approaches  a  constant  as  r 
goes  to  oo.  In  contrast,  the  N3  potential  energy  hypersurface  has  more  complex  structure.  It  has  three 
equivalent  asymptotes  that  look  like  N+N2,  with  the  three  structures  differing  by  the  choice  of  the  lone  N 
atom.  It  also  has  three  local  minima  that  look  like  N-N-N  with  a  bond  angle  of  about  1 15°,  which  differ  in 
the  identity  of  the  central  N  atom,  and  six  local  maxima  which  are  antisymmetric  distortions  of  the  local 
minima  that  look  like  N — N-N  or  N-N — N.  There  are  also  six  global  minima  that  look  like  N«rN-N, 
where  ••  means  van  der  Waals  interaction.  The  van  der  Waals  minima  are  very  shallow.  The  minimum 
energy  pathway  on  the  potential  energy  hypersurface  from  one  N+N2  structure  to  another  involves  passage 
from  the  asymptote  to  the  van  der  Waals  minimum,  to  one  local  maximum,  N — N-N,  to  the  local 
minimum  N-N-N,  to  the  next  local  maximum,  N-N — N,  to  the  van  der  Waals  minimum,  then  on  to  the 
N2+N  asymptote.  The  overall  barrier  height  for  this  process  is  about  3X101  kcal/mol,  which  means  the 
likely  hood  of  N  atom  exchange,  as  this  process  is  called,  is  very  low  except  at  elevated  temperatures 
(several  thousand  degrees).  Finally  the  N3  potential  energy  hypersurface  approaches  a  constant,  high 
energy  as  all  atoms  become  far  apart,  i.e.  forming  3N.  A  depiction  of  the  potential  energy  hypersurface  in 
the  vicinity  of  the  local  maxima  is  shown  in  the  next  figure. 


Figure  5:  N3  PES  for  115°.  Blue  is  low  energy,  green  moderate,  and  red  high  energy. 

For  N+N2  collisions,  there  are  three  possibilities.  The  N  atom  can  simply  bounce  off  the  N2  molecule, 
with  some  energy  transfer  between  translation,  rotation,  and  vibration.  This  is  called  an  inelastic  collision. 
The  second  possibility  is  for  the  trajectory  to  cross  over  from  one  arrangement  valley,  Na+NbNc  in  the 
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figure,  to  another,  say  NaNb+Nc.  This  is  called  an  exchange  reaction,  and  in  contrast  to  inelastic  collisions, 
which  involves  fairly  limited  changes  to  rotational  and  vibrational  energy,  exchange  reaction  products 
tend  to  loose  memory  of  the  initial  arrangement  vibrational  and  rotational  energies  and  lead  to  very  broad 
product  distributions.  The  final  possibility  is  for  the  trajectory  to  climb  to  the  N+N+N  plateau  and 
dissociate  the  N2  molecule. 

How  do  we  describe  these  processes?  This  can  be  done  in  several  ways.  The  most  rigorous  way  is 
to  use  quantum  mechanics.  This  is  quite  expensive,  but  is  feasible  for  energies  below  the  dissociation 
limit.  [3]  Another  way  is  to  use  classical  mechanics  for  some  or  all  of  the  dynamics.  In  our  work,  we  will 
use  the  quasi-classical  trajectory  (QCT)  method.  [8,9]  In  this  method,  we  specify  a  relative  translational 
energy  and  v,j  quantum  numbers,  then  use  old  quantum  mechanics  to  relate  the  quantum  conditions  to  the 
initial  conditions  for  a  classical  trajectory.  The  variables  not  specified  by  Ere\  v,  and  j,  namely  vibrational, 
rotational  phases  and  impact  parameter,  are  integrated  over  using  Monte-Carlo  techniques.  The  final  state 
of  the  system  is  analyzed  using  old  quantum  mechanics.  In  general,  the  final  values  of  the  vibrational  and 
rotational  quantum  numbers,  V  and  will  not  be  integers.  We  turn  them  into  integers  by  rounding  them  to 
the  nearest  integer.  We  thus  obtain  elementary  measure  of  the  outcome  of  a  collision,  the  cross  section 
Opa( Erel),  for  initial  state  a  and  final  state  J3.[ 2]  The  cross  section  has  units  area.  Although  the  calculation 
of  a  single  trajectory  is  not  time  consuming,  the  number  of  trajectories  that  must  be  calculated  is  quite 


Boltzman  distribution 


Ok 


Figure  6:  Energy  distribution  for  averaging  cross  sections. 

large.  For  the  N+N2  calculations  presented  here,  nearly  108  trajectories  were  run. 

Let  us  now  consider  the  issue  of  selecting  the  relative  translational  energy  Erel.  To  compute  a  rate 
coefficient  at  translational  temperature  Tmn\  we  need  to  perform  a  Boltzman  average  of  the  average 
velocity  times  the  cross  section.  In  Figure  6  we  show  the  (un-normalized)  weighting  function  for  the  cross 
section  for  several  temperatures.  At  low  temperature,  the  distribution  is  very  sharply  peaked,  and  the  most 
efficient  way  to  proceed  is  to  sample  energies  around  the  peak  at  each  temperature.  At  higher  temperature, 
the  distribution  becomes  very  broad,  and  it  becomes  more  efficient  to  sample  energies  over  a  wide  range 
and  use  them  for  multiple  temperatures.  Nonetheless,  it  is  not  necessary  to  converge  the  cross  sections  at 
high  energy  as  much  as  for  energies  in  the  vicinity  of  the  peak  in  the  distribution,  for  the  really  high 
energy  cross  sections  have  only  small  impacts  on  the  rate  coefficients.  In  the  present  work,  we  use  a 
combination  of  these  ideas.  We  use  Monte-Carlo  sampling  over  Ere\  with  an  importance  sampling 
function  the  speed  times  the  Boltzman  distribution  at  jtrans= 20000K.  We  then  divide  the  energy  into  64 
intervals,  and  take  the  cross  section  at  the  energy  at  the  middle  of  the  interval  to  be  the  average  cross 
section  we  compute  in  that  energy  interval. 
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In  Figure  7  we  plot  the  computed  dissociation  rate  coefficients  vs.  v  at  Jrran^=20000K.  The  lines 


Figure  7:  Dissociation  rate  coefficients  for  20,000K. 

track  the  v  for  selected  values  of  j.  We  see  that  the  largest  dissociation  rate  coefficient  occurs  at  the  highest 
v  for  each  j ,  with  the  largest  dissociation  rate  being  out  of  j= 0,  v=60.  However,  due  to  the  2j+l  degeneracy 
of  the  rotational  levels, [2]  we  expect  much  larger  populations  for  j  »0,  so  the  j- 0,v=60  level  is  not 
expected  to  play  an  important  role  in  the  dissociation  process.  A  perhaps  more  useful  way  to  understand 
the  dissociation  rate  data  is  to  plot  it  vs.  the  energy  difference  between  the  top  of  the  centrifugal  barrier 
and  the  ro-vibrational  energy  level.  Such  a  plot  is  given  in  the  Figure  8.  Now  we  see  a  much  clearer 
pattern,  with  the  dissociation  rate  showing  a  power  law  relation  with  the  energy  gap. 

2.4  N2+N2 

Now  consider  the  universe  to  consist  of  two  N2  molecules.  Now  the  electronic  structure 
calculations  are  even  more  expensive,  and  calculations  must  be  carried  out  at  thousands  of  geometries. 
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Figure  8:  Dissociation  rates  vs.  dissociation  energy 

Furthermore,  the  topology  of  the  potential  energy  surface  is  quite  complicated.  Below  the  N2+N+N 
asymptote,  there  exist  several  local  minimum.  One  prominent  one  is  the  so-called  tetrahedral  (Td)  N4 
complex.  This  local  minimum  on  the  N4  potential  energy  surface  has  been  proposed  as  a  high  energy 
density  material,  [10]  although  it  has  not  yet  been  synthesized  in  the  laboratory.  The  significance  of  this 

feature  is  if  a  trajectory  samples  this  region,  it  will  loose  memory  of  the  initial  conditions  and  can  lead  to 
rearrangement  scattering,  which  means  unusually  large  vibrational  and  rotational  energy  transfer.  There 
are  also  other  local  minima,  such  as  the  rectangle  form  of  N4,  of  which  there  are  a  number  of  equivalent 
forms.  This  all  makes  the  determination  of  a  reliable  analytic  representation  very  difficult.  Our  best 
strategy  so  far  is  illustrated  in  Figure  9.  There  we  show  results  of  electronic  structure  calculations  for  the 
cross  or  X  orientation  as  a  function  of  one  N2  bond  length  for  two  different  values  of  the  N2-N2  distance. 
The  symbols  are  color  coded  so  that  green  means  the  N2  bond  length  is  smaller  than  the  N2-N2  distance,  so 
the  system  is  more  like  N2+N2,  while  red  means  larger  N2  bond  lengths,  so  the  system  is  more  like 
N2+N+N.  The  dark  purple  curve  is  the  result  of  diagonalizing  a  4x4  matrix,  with  diagonals  analytic 
representations  of  potential  energies  of  systems  that  look  like  N2+N2,  N2+N+N,  Td  N4,  and  rectangular  N4, 
and  off  diagonals  being  constants.  For  the  small  N2-N2  distance,  we  also  show  the  Td  N4  (light  blue)  and 
N2+N+N  (purple)  diagonals.  Although  this  analytic  representation  is  still  only  qualitatively  correct,  one 
can  clear  see  how  the  different  topological  features  of  the  potential  come  into  play. 

For  nuclear  dynamics,  we  can  follow  the  same  procedure  as  for  N+N2.  Now,  however,  quantum 
mechanical  dynamics  calculations  will  be  possible  only  for  quite  low  energy.  In  contrast,  the  QCT  method 
can  be  used  without  much  more  effort  per  trajectory  compared  to  N+N2.  Once  we  have  finalized  out 
N2+N2  potential,  we  will  be  carrying  out  trajectories  for  N2+N2  on  a  similar  scale  as  for  N+N2. 
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Figure  9:  Cuts  of  N4  potential  in  X  orientation. 


2.4  Gaseous  N,  N2 

We  now  transition  from  the  atomic  level  to  a  macroscopic  level.  Once  again  we  can  take 
advantage  of  the  properties  of  nature  to  simplify  our  simulations.  As  we  saw,  as  we  went  from  two,  to 
three,  to  four  atoms,  our  calculations  got  much  more  difficult.  This  will  continue  as  we  go  from  four  to  16 
to  32  ...  But  how  far  to  we  have  to  go?  The  unit  of  measure  of  the  amount  of  a  substance  at  the 
macroscopic  level  is  Avogadro's  number,  NA=6.02xl023.  This  is  the  number  of  water  molecules  in  a  little 
more  that  one  Tablespoon.  This  same  amount  of  water,  if  allowed  to  evaporate  at  room  temperature  and 
pressure,  would  occupy  about  22.4  liters.  Avogadro's  number  is  huge,  and  it  will  not  be  possible  to  treat 
this  number  of  individual  molecules.  Fortunately,  the  fact  that  NA  is  so  large  also  means  that  statistical 
treatments  will  be  very  accurate.  [2]  This  means  that  the  nature  of  our  questions  change,  that  is,  rather  than 
asking  what  the  outcome  of  a  collision  is,  we  ask  how  level  populations  change  with  time.  The  equation 
we  will  use  to  describe  this  is  the  master  equation. 

In  the  master  equation  formulism,  we  assume  that  translational  degrees  of  freedom  equilibrate 
instantaneously  to  a  distribution  governed  by  the  temperature  Ttrans.  See  [11]  for  an  example  of  previous 
work  by  Schwenke  on  the  H,  H2  system.  The  translational  temperature  can  be  fixed,  or  driven  by  coupling 
to  flow  equations.  The  solution  then  describes  the  time  dependence  of  the  concentrations  for  the  processes 


N20»^2N 

N20»  +  NoN2C/',v')  +  N 

N20»  +  N2(/,v)  N2(/,v')  +  N2(7,v') 

N2(/,v)  +  No3N  (1> 

N20»  +  N2(7,v)<h>N2(/>')  +  2N 
N2(/,v)  +  N2(/,v)^4N 

for  external  temperature  Ttrans ,  where  j  is  the  rotational  quantum  number  and  v  is  the  vibrational  quantum 
number.  The  master  equation  can  be  written  in  many  forms,  and  the  one  we  choose  for  Eq.  (1)  is 
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where  njv  is  the  concentration  of  N 2(/,v),  nN  is  the  concentration  of  N  atoms,  the  sums  are  over  all  bound 
and  all  long  lived  metastable  levels  of  N2,  and  the  k  are  pseudo  first  or  second  order  rate  coefficients 
defined  as 
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where  5nm  is  unity  if  n=m  and  zero  otherwise, 


kET 


is  an  energy  transfer  rate  coefficient, 


kD  is  an 


Figure  10:  populations  for  H,  H2  from  master  equation. 

dissociation  rate  coefficient,  kR  is  an  recombination  rate  coefficient,  kPD  is  a  predissociation  rate 
coefficient,  kTR  is  a  tunnelling  recombination  rate  coefficient,  kDD  is  a  double  dissociation  rate  coefficient, 
and  kDR  is  a  double  recombination  rate  coefficient.  The  predissociation  and  tunnelling  recombination  rate 
coefficients  are  nonzero  only  for  metastable  jv  levels.  Since  N2  has  9390  jv  levels,  this  means  for  the  N,  N2 
system  there  will  9391  coupled  ordinary  differential  equations,  thus  the  solution  of  the  master  equation 
will  be  very  arduous. 

For  an  example  of  what  can  be  learned  from  a  master  equation  study  coupled  with  accurate  state  - 
to-state  rate  coefficients,  readers  are  referred  to  Ref.  11.  Here  we  studied  the  H,  H2  system.  This  system 
shares  a  great  deal  with  the  N,  N2  system,  except  the  fact  that  H2  has  only  348  jv  levels  makes  the  solution 
of  the  master  equation  much  easier.  In  Figure  10  we  show  populations  as  a  function  of  time  for  a  2000K 
temperature  jump.  This  figure  is  a  distillation  of  a  vast  amount  of  data,  yet  we  see  some  very  interesting 
aspects  of  the  dissociation  of  H2.  First  we  note  that  the  average  value  of  j  and  v  equilibrate  at  about  the 
same  rate,  reaching  their  final  values  by  a  few  psec.  In  contrast,  dissociation  does  not  begin  to  occur  until 
about  10  psec,  with  vibrational  and  rotation  in  equilibrium.  Now  the  different  vibrational  rotational 
spacings,  potential  energy  surfaces,  and  temperatures  for  the  N,N2  system  will  be  different  from  the  H,  H2 
system,  so  this  decoupling  of  time  scales  may  not  occur  in  our  system,  but  this  remains  a  very  tantalizing 
result. 


A  serious  difficulty  addressed  in  Ref.[ll]  was  how  to  compare  the  master  equation  simulations  to 
experiment.  This  is  because,  e.g .,  experimental  measurements  report  the  rate  kD  in  the  expression 


dn v- 


dt 
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and  the  master  equation  can  not  be  analytically  reduced  to  this  form.  Therefore  kD  is  known  as  the 


David  W.  Schwenke:  H 


FTG.  8.  Comparison  of  the  current  results  for  k 3bdy  Hj  (□),  to  previous 
(Ref.  1 )  ORT  results  (O)  and  the  experimental  recommendation  of  Ref.  3 
(solid  line).  The  dotted  line  demarcates  the  experimental  error  limit. 


Figure  11:  Comparison  between  calculated  and  experimental  recombination  rate  coefficients 
for  H2. 

phenomenological  rate  coefficient  because  although  the  rate  law  of  Eq.  (8)  is  observed,  it  can  not  be 
derived  from  more  accurate  theory.  To  compare  to  experiment,  we  considered  methods  based  on  three 
broad  ideas:  (i)  data  fitting,  (ii)  equilibrium  one  way  flux,  and  (iii)  eigenvalue  analysis. 
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In  the  data  fitting  approach,  we  closely  mimicked  experiment  by  determining  the  time  dependence 
of  the  N  and  aggregate  N2  concentrations,  and  then  determining  the  best  kD  that  reproduced  the  time 
profile  when  we  used  Eq.  (8).  In  Figure  1 1  we  show  the  comparison  with  experiment  that  we  obtained. 

The  results  were  very  good.  It  should  be  noted  that  the  apparent  deterioration  at  the  highest  temperature  is 
just  simply  a  reflection  on  the  uncertainty  in  determining  the  rate  coefficient  since  at  these  temperatures, 
the  equilibrium  faction  of  H2  is  very  small. 

The  next  scheme  we  used  to  extract  the  phenomenological  rate  coefficient  was  the  equilibrium  one¬ 
way  flux  method.  This  was  not  nearly  as  reliable  a  method,  being  only  capable  of  order  of  magnitude 
predictions  of  the  recombination  rate  coefficient.  Nevertheless  this  will  give  valuable  insight  into  the 
dissociation  process.  In  this  method,  we  evaluate  the  rhs  of  Eq.  (2)  using  equilibrium  concentrations  and 
only  include  terms  leading  to  dissociation.  We  show  the  contribution  of  the  various  terms  as  a  function  of 
v  in  Figure  13  and  as  a  function  of  j  in  Figure  12.  We  see  that  the  dissociation  rate  is  dominated  by  low  v 
and  high  j. 


Figure  12:  Equilibrium  dissociation  flux  at  20,000K. 
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Figure  13:  Equilibrium  dissociation  flux  at  20,000K. 

We  next  turn  to  the  eigenvalue  analysis  method.  For  H2,  it  proved  to  be  a  very  accurate  and  efficient 
method.  In  this  method,  we  observe  that  we  can  write  the  master  equation  in  the  form 

dn  T  r 

—  =  C(n)n  =  C{nf)n  + ...  (9) 

where  bold  means  matrix  and  the  superscript /means  final.  That  is,  we  linearize  the  master  equation.  Then 
the  solution  can  be  found  in  terms  of  the  eigenvalues  and  eigenvectors  of  the  matrix  C: 

"«(0  =  ni  +YaVaPCp  eXP  V  (10) 

P 

Then  the  phenomenological  rate  coefficient  is  computed  from  the  first  nonzero  eigenvalue.  The  problem 
with  this  procedure  for  the  N2  system  is  that  C  is  a  9391x9391  non-symmetric  matrix,  and  good  canned 
routines  to  determine  the  lowest  eigenvalues  and  eigenvectors  of  a  matrix  of  this  size  do  not  exist.  But  this 
can  be  overcome  by  the  following  technique:  if  we  multiply  C  from  the  right  by  the  diagonal  matrix  made 
up  of  the  square  roots  of  the  nfa  and  then  multiply  the  result  from  the  left  with  the  inverse  of  this  diagonal 
matrix,  then  we  obtain  a  new  matrix  which  is  made  up  of  a  9390x9390  symmetric  part,  and  then  one  row 
and  one  column  which  are  not  symmetric.  We  can  then  diagonalize  the  symmetric  part  efficiently, 
determining  only  the  eigenvalues  of  interest,  and  then  use  the  non-symmetric  Jacobi  method  to  clean  up 
the  rest  to  yield  the  eigenvalues  of  interest.  Thus  we  expect  this  to  be  a  very  promising  method  for  treating 
the  full  system. 


3.0  MODELING 

Finally  let  us  ask  the  question  of  what  is  it  that  computational  fluid  dynamics  actually  needs  for 
accurate  simulations.  Clearly  following  the  populations  of  all  9390  vj  levels  is  a  vast  overkill  and  also 
prohibitatively  expensive.  What  is  actually  needed  are  the  N,  and  N2  number  densities,  to  go  into  the  mass 
conservation  equations,  the  average  internal  energies,  to  go  into  the  energy  conservation  equations,  and 
finally  the  internal  state  distribution,  to  go  into  the  radiation  transfer  equations.  The  eigenvalue  analysis 
method  gives  all  this  information,  except  the  eigenvalues  and  eigenvectors  depend  on  Tmn\  the  total 
density  and  the  particular  mixture.  Thus  a  more  flexible  method  is  required. 

Let  us  suppose  that  all  the  vj  state  number  densities  are  not  independent,  i.e. 
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(11) 


na(t)  =  x1(t)n°a,  a  =  l,...,N1 

na(t)  =  x2(t)n°a,  a=Nl+l,...,N2 
etc.  Then  can  we  determine  how  the  v  depend  on  time?  This  technique  is  known  as  lumping,  and  can  be 
written  mathematically  as  follows:  let 

x  =  M n  (12) 

where  M  is  some  invertible  9391x9391  matrix.  Then  the  governing  equation  for  the  v  is 

dx  r 
—  =  Kx 
dt 

with 


(13) 


K=MCM'  (14) 

So  far  this  has  just  introduced  a  large  amount  of  extra  work,  but  now  make  the  approximation  that  xt= 0  for 
all  time  for  i>m,  m«9391.  This  converts  the  master  equation  into  one  having  only  m  coupled  equations. 
The  advantage  of  this  method  is  by  performing  a  significant  amount  of  work  up  front  to  determine  the  rate 
equations  for  the  lumps,  we  can  recover  accurate  time  dependent  populations  from  a  much  smaller  set  of 
equations. 


What  are  some  ways  to  choose  the  lumps?  One  choice  is  to  have  m  run  0  to  60  and  represent  the 
vibrational  levels  with  the  j  populations  be  Boltzman  at  some  temperature,  say  Trans.  This  is  a  perfectly 
well  defined  and  reasonable  proposition,  provided  the  j  levels  are  those  for  the  particular  v.  Thus  we  can 
eliminate  the  j  levels,  as  we  would  like  to  do,  but  note  that  this  is  not  the  same  is  decoupling  vibration  and 
rotation.  Some  other  lumping  schemes  are  to  sort  jv  levels  by  energy,  and  include  in  a  lump  all  levels  with 
similar  energies.  This  has  the  advantage  that  the  populations  within  the  lump  will  be  primarily  determined 
by  degeneracy  factors,  rather  than  an  external  temperature.  Another  idea  would  be  to  carry  out  the 
eigenvalue  analysis  for  a  series  of  likely  situation,  and  then  use  the  eigenvectors  to  define  the  lumps.  And 
of  course,  other  ideas  are  possible. 

4.0  CONCLUSIONS 

I  have  outlined  the  procedure  involved  from  going  from  a  first  principle  treatment  of  molecular 
dissociation  to  a  model  for  use  in  CFD  calculations.  There  are  a  multitude  of  complex  problems  that  must 
be  solved  before  a  viable  model  is  developed.  Our  goal  is  to  use  the  first  principle  data  to  determine  a 
hierarchal  set  of  models  that  can  be  used  in  CFD  calculations.  These  models  will  have  wide  applicability, 
and  systematically  trade  off  between  accuracy  and  expense. 
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